Propaganda Wars
Amanda Franklin-Ryan
2023-04-25
Introduction
Atrocity Propaganda
“History stopped in 1936- after that, there was only propaganda” - George Orwell
“We are children no more”, 15-year-old ‘Nayirah’” tells the Congressional Caucus, blinking back tears, her voice cracking as she describes Iraqi soldiers storming the hospital where she volunteered, stealing incubators and throwing babies onto the floor to die. It’s 1991, Vietnam’s legacy still looms large and the American public have no desire to be embroiled in yet another distant war, but Nayirah’s testimony sways them with the US invading Iraq 3-months later. Except it wasn’t true, any of it. Nayirah al-Sabah, the daughter of the Kuwaiti ambassador had never set foot in the hospital, instead she was coached to perfection for her Oscar winning performance by PR firm Hill and Knowlton who were directly employed by the Kuwaiti Government.(Abuhamad 1994) (“How False Testimony and a Massive U.S. Propaganda Machine Bolstered George H.W. Bush’s War on Iraq,” n.d.)
Atrocity propaganda is nothing new, with lurid tales of fabricated massacres stretching back to the Crusades, but now storytelling is easier than ever before; all you need is a mobile and internet connection. There are now more hours of Syrian conflict footage online, than time elapsed since the war started.Almost five million people are currently trapped in Khartoum as Sudanese generals Abdel Fattah al-Burhan and Mohamed Hamdan Daglo vie for power and Twitter is awash with fake footage taken from conflicts in Yemen, Libya and even a computer game. With most Twitter accounts geolocation data switched off, it’s game on for opportunists worldwide to create and share their own realities. Even with image metadata and reverse image search, it’s almost impossible to ascertain whether shared images are accurate, misrepresentative, staged or concocted by government backed bots. DALL.E2 will only make things worse, with the potential power to conjure up bespoke atrocities in every style from gritty World War II grayscale to heartrending full colour closeup at the click of a mouse.
Yet in a world where fiction seems higher definition than reality, satellite imagery could offer a glimmer of hope. In March 2022 Ukrainians reported the rape, murder and torture of hundreds of civilians in Bucha, with bodies strewn across the streets. Unsurprisingly, the Russian Defence Ministry immediately dismissed the allegations as, “another production of the Kyiv regime for the Western media”, blaming Ukrainian soldiers themselves for the killings. However Maxar satellite data corroboratedthe Ukrainian narrative illustrating the bodies had been present long before Ukrainian forces entered the area.
While high resolution Maxar data is only available commercially or for selected humanitarian disasters, medium resolution Sentinel-1 and Sentinel 2 satellite imagery is freely accessible online. Evidently it’s no magic bullet. evidence of newly dug mass graves or bombed buildings, proves only the incident is not a recycled copy-paste of a past atrocity and in isolation can rarely identify the true perpetrator, however it is already being used in conflict mapping. Depending on the location and nature of the violence, burning oil wells; abandoning agricultural land, (pictures here) construction of refugee camps and building damage have been used to identify and monitor armed conflict. Using methods to map meteor craters on the moon, a 2020 studyused satellite imagery combined with machine learning to identify unexploded munitions from the Vietnam War embedded in overgrown bomb craters in rural Cambodia. Detecting chemical and biological warfare with satellite imagery is less promising due to the low concentrations of the agents in the atmosphere, however vegetation changes could potentially be use as a proxy for contamination with radioactivity and journalists are already using satellite pictures to investigate North Korea’s chemical weapon capabilities.
Raqqa: The Perfect Perpetratorless Crime
Created in 2019, Amnesty International’s Strike Tracker used satellite imagery coupled with crowd sourced image labelling to create an accurate timeline for the 2017 US led Coalition bombing in the Syrian city of Raqqa. Attempting to destroy Islamic State’s “Caliphate”, Coalition forces fired over 30,000 missiles in a four-month bombardment, destroying 80% of the city and killing and injuring thousands of civilians. Pinning down responsibility was challenging as the Coalition accused ISIS, Russia blamed the US and the US denied all responsibility ruling less than one percent of civilian casualty reports “credible”.(April 25 2019 and P.m, n.d.)
Amnesty International enlisted thousands of digital activists to analyse a series of satellite photos of the same building to identify the exact date it was destroyed. After following a short interactive tutorial volunteers were presented with a new series of images to label. Combining the crowd sourced strike timeline with witness testimony and interviews, Amnesty forced the Coalition to admit responsibility for 159 civilian deaths.(“IS Conflict: Coalition Strikes on Raqqa ’Killed 1,600 Civilians’” 2019)
This project aims to collect the dataset necessary to replicate Amnesty’s Strike Tracker Project. However instead of using high resolution human labelled satellite images, medium resolution Sentinel-1 and Sentinel-2 images will be downloaded and in a future seminar paper (hopefully) used to investigate the efficacy of using machine learning methods to automatically detect building damage and to generate a destruction timeline using a simplified version of the methodology described in x. More specifically the project will:
Collect and visualise data from UNOSAT illustrating the position of the coalition strikes (for the purposes of the future machine learning model this represents the ground truth)
Identify and download suitable Sentinel-2 satellite images from the SentinelHub API, an R function will also be written to simplify this process
Create a GIF animation with the downloaded satellite data illustrating the destruction of Raqqa
Position of Coalition Strikes: UNOSAT Data
Since 2001 the United Nations Satellite Centre (UNOSAT) has published and analysed geospatial data on climate change, natural disaster and conflict including this 2017 dataset on the buildings destroyed by the US-Coalition bombing of Raqqa The shape file contains data taken from five different sets of satellite images from (2013-2015) and includes information on the level of damage (moderate, severe, destroyed); position of the damaged building (latitude, longitude, neighbourhood) and the confidence level in the degree of destruction (Very high, medium, uncertain). It also contains a SiteID column with further information on the nature of the destroyed building (church, university, stadium etc) However the vast majority of buildings are simply labelled “Building (General / Default)”. Building destruction has been manually labelled using imagery from the Pleiades, Worldview-2, Worldview-1 and Worldview-3 satellites and has not been verified on the ground.
library(rio)
library(data.table)
summary <- rio::import("../Visualisations/head_original_raqqa_data.rds")
summary_table <- data.table(summary)
summary_table## SiteID SensDt SensId Confidence
## 1: Building (General / Default) 2013/10/22 Worldview-1 Very High
## 2: Building (General / Default) 2013/10/22 Worldview-1 Very High
## 3: Building (General / Default) 2013/10/22 Worldview-1 Very High
## 4: Building (General / Default) 2013/10/22 Worldview-1 Medium
## 5: Building (General / Default) 2013/10/22 Worldview-1 Very High
## 6: Building (General / Default) 2013/10/22 Worldview-1 Very High
## 7: Building (General / Default) 2013/10/22 Worldview-1 Very High
## 8: Building (General / Default) 2013/10/22 Worldview-1 Very High
## 9: Building (General / Default) 2013/10/22 Worldview-1 Very High
## 10: Building (General / Default) 2013/10/22 Worldview-1 Medium
## DaSitCl SensorDate SensID2 ConfID2 DaSitCl2
## 1: Destroyed 2014/02/12 Worldview-2 Very High Destroyed
## 2: Destroyed 2014/02/12 Worldview-2 Very High Destroyed
## 3: Destroyed 2014/02/12 Worldview-2 Very High Destroyed
## 4: Moderate Damage 2014/02/12 Worldview-2 Medium Moderate Damage
## 5: Destroyed 2014/02/12 Worldview-2 Very High Destroyed
## 6: Destroyed 2014/02/12 Worldview-2 Very High Destroyed
## 7: Destroyed 2014/02/12 Worldview-2 Very High Destroyed
## 8: Severe Damage 2014/02/12 Worldview-2 Medium No Visible Damage
## 9: Severe Damage 2014/02/12 Worldview-2 Medium No Visible Damage
## 10: Severe Damage 2014/02/12 Worldview-2 Medium Severe Damage
## Damst2 SensDt3 SensID3 ConfId3 DaSitCl3
## 1: No change 2015/05/29 Pleiades Very High Destroyed
## 2: No change 2015/05/29 Pleiades Very High Destroyed
## 3: No change 2015/05/29 Pleiades Very High Destroyed
## 4: No change 2015/05/29 Pleiades Medium No Visible Damage
## 5: No change 2015/05/29 Pleiades Very High Destroyed
## 6: No change 2015/05/29 Pleiades Very High Destroyed
## 7: No change 2015/05/29 Pleiades Very High Destroyed
## 8: Decrease - damage 2015/05/29 Pleiades Medium No Visible Damage
## 9: Decrease - damage 2015/05/29 Pleiades Medium No Visible Damage
## 10: No change 2015/05/29 Pleiades Medium Severe Damage
## Damst3 SensDt4 SensID4 ConfId4 DaSitCl4
## 1: No change 2017/02/03 Worldview-2 Very High Destroyed
## 2: No change 2017/02/03 Worldview-2 Very High Destroyed
## 3: No change 2017/02/03 Worldview-2 Very High Destroyed
## 4: Decrease - damage 2017/02/03 Worldview-2 Medium No Visible Damage
## 5: No change 2017/02/03 Worldview-2 Very High Destroyed
## 6: No change 2017/02/03 Worldview-2 Very High Destroyed
## 7: No change 2017/02/03 Worldview-2 Very High Destroyed
## 8: No change 2017/02/03 Worldview-2 Medium No Visible Damage
## 9: No change 2017/02/03 Worldview-2 Medium No Visible Damage
## 10: No change 2017/02/03 Worldview-2 Medium Severe Damage
## Damst4 SensDt5 SensId5 ConfId5 DaSitCl5 Damst5
## 1: No change 2017/10/21 WorldView-3 Very High Destroyed No change
## 2: No change 2017/10/21 WorldView-3 Very High Destroyed No change
## 3: No change 2017/10/21 WorldView-3 Very High Destroyed No change
## 4: No change 2017/10/21 WorldView-3 Medium No Visible Damage No change
## 5: No change 2017/10/21 WorldView-3 Very High Destroyed No change
## 6: No change 2017/10/21 WorldView-3 Very High Destroyed No change
## 7: No change 2017/10/21 WorldView-3 Very High Destroyed No change
## 8: No change 2017/10/21 WorldView-3 Medium No Visible Damage No change
## 9: No change 2017/10/21 WorldView-3 Medium No Visible Damage No change
## 10: No change 2017/10/21 WorldView-3 Medium Severe Damage No change
## GrDaCl FieldVl Notes Name Neigh Code
## 1: Damaged Buildings Not yet field validated <NA> Raqqa Andalus CE20130604SYR
## 2: Damaged Buildings Not yet field validated <NA> Raqqa Andalus CE20130604SYR
## 3: Damaged Buildings Not yet field validated <NA> Raqqa Andalus CE20130604SYR
## 4: Damaged Buildings Not yet field validated <NA> Raqqa Andalus CE20130604SYR
## 5: Damaged Buildings Not yet field validated <NA> Raqqa Andalus CE20130604SYR
## 6: Damaged Buildings Not yet field validated <NA> Raqqa Andalus CE20130604SYR
## 7: Damaged Buildings Not yet field validated <NA> Raqqa Andalus CE20130604SYR
## 8: Damaged Buildings Not yet field validated <NA> Raqqa Andalus CE20130604SYR
## 9: Damaged Buildings Not yet field validated <NA> Raqqa Andalus CE20130604SYR
## 10: Damaged Buildings Not yet field validated <NA> Raqqa Andalus CE20130604SYRData has only been collected once a year from 2013-2016 and twice in 2017 so it is impossible to know the exact date a building was destroyed,. all that can be established is the date of the first satellite image in which it was damaged. 13155 buildings were damaged over the entire timeframe, a few hundred in ad hoc airstrikes from 2013-2016, but over 10,000 in 2017. Evidently in the earlier years, buildings would have been damaged and repaired (possibly multiple times) but for the purposes of this analysis, a building can only be labelled “destroyed”,“severely damaged” etc. once at the first instance of damage. The dataframe was reshaped from a geospatial frame to dataframe to reflect this using the code below. In this case a spatial dataframe wasn’t required as no polygon data was included.
library(rgdal) # loading spacial data
library(sp)
library(raster)
library(dplyr)
library(data.table)
# Load and view dataset
directory <- "UNOSAT Data/RaqqaDamage"
raqqa_data <- rgdal::readOGR(directory,"Damage_Sites_Raqqa_CDA")
head_raqqa_data <- head(raqqa_data,10)
unique(raqqa_data$ConfId5)
head_raqqa_data_table <- as.data.table(head_raqqa_data)
formatted_head_raqqa <- formattable(head_raqqa_data_table,
align=c("l","l"))
saveRDS(formatted_head_raqqa,"Visualisations/head_original_raqqa_data.rds")
raqqa_data@data$SiteID
###---- 1 Identify key characteristics of spacial data
str(raqqa_data)
# Data from 5 sets of satellite images from 2013-2017
# For each set of satellite data contains information on level and destruction and level of confidence.
# Estimates of destruction made entirely from satellite imagery and not verified on the ground
#Data projection
raqqa_data@proj4string
# No projection just latitude and longitude WG584
# Coordinates of area covered by dataset
raqqa_data@bbox
#Coordinates of damaged buildings
raqqa_data@coords
#plot raqqa data
plot(raqqa_data)
###--- Reshape Data ---
# Reshape data to find year when buildings were identified as "Destroyed", "Moderate Damage" or "Severe Damage"
# Buildings can remain destroyed over time or be rebuilt, it's important not to double count
# For the purposes of this analysis the same building can only be classified as damaged/destroyed once
# even if in reality it has been bombed and rebuilt several times, this could lead to under counting destruction
raqqa_dataframe <- as.data.frame(raqqa_data)
# Remove duplicate rows
raqqa_dataframe <-raqqa_dataframe %>%
distinct()
#Identify first date building destruction detected
destroyed_dates <- raqqa_dataframe %>%
mutate(damage_year=case_when(DaSitCl=="Destroyed"~"2013/10/22",
DaSitCl2=="Destroyed"~"2014/02/12",
DaSitCl3=="Destroyed"~"2015/05/29",
DaSitCl4=="Destroyed"~"2017/02/03",
DaSitCl5=="Destroyed"~"2017/10/21"))
#Add destroyed column
destroyed_dates <- destroyed_dates %>%
mutate(destruction_type=if_else(!is.na(damage_year),"Destroyed",damage_year))
# Extract destroyed entries into single dataframe
destroyed_dataframe <- destroyed_dates %>%
filter(destruction_type=="Destroyed")
# Remove rows from data frame and use same method to identify dates buildings had "Moderate Damage"
severe_damaged_dates <- destroyed_dates %>%
filter(is.na(destruction_type)) %>% # Can't use != as doesn't handle NAs well
mutate(damage_year=case_when(DaSitCl=="Severe Damage"~"2013/10/22",
DaSitCl2=="Severe Damage"~"2014/02/12",
DaSitCl3=="Severe Damage"~"2015/05/29",
DaSitCl4=="Severe Damage"~"2017/02/03",
DaSitCl5=="Severe Damage"~"2017/10/21"))
#Add "Severe damage" column
severe_damaged_dates <- severe_damaged_dates %>%
mutate(destruction_type=if_else(!is.na(damage_year),"Severe Damage",damage_year))
# Extract "Severe damage" into single dataframe
severe_damage_dataframe <- severe_damaged_dates %>%
filter(destruction_type=="Severe Damage")
# Remove rows from data frame and use same method to identify dates buildings had "Moderate Damage"
moderate_damaged_dates <- severe_damaged_dates %>%
filter(is.na(destruction_type)) %>% # Can't use != as doesn't handle NAs well
mutate(damage_year=case_when(DaSitCl=="Moderate Damage"~"2013/10/22",
DaSitCl2=="Moderate Damage"~"2014/02/12",
DaSitCl3=="Moderate Damage"~"2015/05/29",
DaSitCl4=="Moderate Damage"~"2017/02/03",
DaSitCl5=="Moderate Damage"~"2017/10/21",
SiteID=="Road"~"2017/10/21",
SiteID=="Field"~"2017/10/21"))
#Add "Moderate damage" column
moderate_damaged_dates <- moderate_damaged_dates %>%
mutate(destruction_type=if_else(!is.na(damage_year),"Moderate Damage",damage_year))
# Extract Moderate" entries into single dataframe
moderate_damage_dataframe <- moderate_damaged_dates %>%
filter(destruction_type=="Moderate Damage")
# Combine destroyed, severe and moderate dataframes into single dataframe
raqqa_dataframe_reshaped <- rbind(severe_damage_dataframe, moderate_damage_dataframe, destroyed_dataframe)
saveRDS(raqqa_dataframe_reshaped,"Datasets/Reshaped UNOSAT data.rds")For this analysis, I decided to focus exclusively on the coalition airstrikes of 2017 and was interested in visualising the type and extent of the damage, I experimented with bar charts and summaries of the number of each building destroyed (see code below), however concluded that an interactive map would be the best way to present the information using R’s leaflet package.
library(dplyr)
library(ggplot2)
library(lubridate)
library(data.table)
library(formattable)
library(sen2r)
#Load UNOSAT data in dataframe form
raqqa_data <- readRDS("Datasets/Reshaped UNOSAT data.rds")
raqqa_data_summarised <- raqqa_data %>%
select(SiteID,Neigh,damage_year,destruction_type,coords.x1,coords.x2) %>%
mutate(damage_year=lubridate::ymd(damage_year))
### --- Exploratory Data Analysis---
#Total number of buildings destroyed
nrow(raqqa_data_summarised)
# 13155 buildings destroyed over 5 year period
### 1. When were most buildings destroyed or damaged?
no_building_destroyed <- raqqa_data_summarised %>%
group_by(damage_year, destruction_type) %>%
summarise(n())
#Unsurprisingly most damage occurred during coalition bombing in October 2017
# Over 100000 buildings damaged during this time
# Bar chart illustrating severity of damage by date
year_labels <- c("October 2013","February 2014","May 2015","February 2017", "October 2017")
damage_plot <- ggplot(data=no_building_destroyed, aes(x=damage_year, y=`n()`, fill=destruction_type)) +
geom_bar(stat="identity", position=position_dodge())+
theme_minimal()+
ylab("Number of buildings")+
xlab("")+
scale_x_discrete( labels= year_labels)+
ggtitle("Severity of damage to buildings in Raqqa (2013-2017)")
### 2. What kind of buildings were destroyed?
# Calculate type of buildings destroyed by year
buildings_destroyed <- raqqa_data_summarised %>%
group_by(damage_year,SiteID) %>%
summarise(n()) %>%
arrange(n()) %>%
ungroup()
# Focusing exclusively on the coalition bombing of October 2017
coalition_bombing <- buildings_destroyed %>%
filter(damage_year== "2017-10-21") %>%
arrange(desc(`n()`))
names(coalition_bombing) <- c('Year','Building type','Number')
coalition_bombing <- coalition_bombing %>%
select(`Building type`, Number)
coalition_bombing_data_table <- as.data.table(coalition_bombing)
formatted_table <- formattable(coalition_bombing,
align=c("l","l"))
#3. Which neighbourhoods were targeted?
neighbourhoods_targeted <- raqqa_data_summarised %>%
group_by(damage_year,Neigh) %>%
summarise(sum(total=n())) %>%
arrange(sum(total=n())) %>%
ungroup()
# What are the coordinates for the bounding box around Raqqa (needed for download of Sentinel data))
max_x1 <- max(raqqa_data_summarised$coords.x1)
min_x1 <- min(raqqa_data_summarised$coords.x1)
max_x2 <- max(raqqa_data_summarised$coords.x2)
min_x2 <- min(raqqa_data_summarised$coords.x2)
# What are the coordinates of the "centre" of the city?
centre_x1 <- mean(raqqa_data_summarised$coords.x1)
centre_x2 <- mean(raqqa_data_summarised$coords.x2)library(leaflet)
library(dplyr)
library(rio)
library(leaflet.extras) # for extra features like adding search box and reset map button
library(htmltools) #for escaping html
###--- Create interactive leaflet map of Raqqa
###--- Import UNOSAT data and focus on October 2017 (Coalition bombing)
unosat_data <- rio::import("../Datasets/Reshaped UNOSAT data.rds")
coalition_strikes_2017 <- unosat_data %>%
filter(damage_year=="2017/10/21")
###--- Create plot illustrating distribution of destroyed, moderately damaged and destroyed buildings---
# Create custom destruction colour palette
building_palette <- leaflet::colorFactor(palette=c("#00876c","#003f5c", "#85b96f", "#bdcf75","#f6bc63","#f19452","#e66a4d","#d43d51"),
levels=c(" General Building","Field","Road","School/University","Industrial Facility","Mosque","Government Building","Hospital"))
###--- Create plot with groups to make it look less cluttered---
# Dataframe with only generic buildings buildings
building_only_dataframe <- coalition_strikes_2017 %>%
filter(SiteID=="Building (General / Default)")
grouped_building_plot <- Building_only_dataframe %>%
leaflet()%>%
addProviderTiles("OpenStreetMap") %>%
setView(lng=39.00897, lat= 35.95394, zoom=14) %>%
setMaxBounds(lng1 = 39.05947,
lat1 = 35.9745,
lng2 = 38.93641,
lat2 = 35.92829) %>%
clearMarkers() %>%
addCircleMarkers(lng=building_only_dataframe$coords.x1,
lat=building_only_dataframe$coords.x2,
radius = 2,
stroke=FALSE,
fillOpacity= 0.8,
color = "#4fa16e",
label=~destruction_type,
group="General Building")
# Dataframe with only Fields
field_only_dataframe <- coalition_strikes_2017 %>%
filter(SiteID=="Field")
#Create field layer
grouped_building_plot <- grouped_building_plot %>%
addCircleMarkers(lng=field_only_dataframe$coords.x1,
lat=field_only_dataframe$coords.x2,
radius = 2,
stroke=FALSE,
fillOpacity= 1,
color = "#003f5c",
label=~destruction_type,
group="Field") %>%
addLayersControl(overlayGroups = c("General Building", "Field"))
# Dataframe with only roads
road_only_dataframe <- coalition_strikes_2017 %>%
filter(SiteID=="Road")
#Create road only layer
grouped_building_plot <- grouped_building_plot %>%
addCircleMarkers(lng=road_only_dataframe$coords.x1,
lat=road_only_dataframe$coords.x2,
radius = 2,
stroke=FALSE,
fillOpacity= 1,
color = "#85b96f",
label=~destruction_type,
group="Road") %>%
addLayersControl(overlayGroups = c("General Building", "Field", "Road"))
# Dataframe with only schools/universities
school_only_dataframe <- coalition_strikes_2017 %>%
filter(SiteID=="School / University")
#Create school only layer
grouped_building_plot <- grouped_building_plot %>%
addCircleMarkers(lng=school_only_dataframe$coords.x1,
lat=school_only_dataframe$coords.x2,
radius = 2,
stroke=FALSE,
fillOpacity= 1,
color = "#bdcf75",
label=~destruction_type,
group="School/University") %>%
addLayersControl(overlayGroups = c("General Building", "Field", "Road", "School/University"))
# Dataframe with only Industrial facilities
industrial_only_dataframe <- coalition_strikes_2017 %>%
filter(SiteID=="Industrial Facility")
#Create industrial only layer
grouped_building_plot <- grouped_building_plot %>%
addCircleMarkers(lng=industrial_only_dataframe$coords.x1,
lat=industrial_only_dataframe$coords.x2,
radius = 2,
stroke=FALSE,
fillOpacity= 1,
color = "#f6bc63",
label=~destruction_type,
group="Industrial Facility") %>%
addLayersControl(overlayGroups = c("General Building", "Field", "Road", "School/University", "Industrial Facility"))
# Dataframe with only mosques
mosque_only_dataframe <- coalition_strikes_2017 %>%
filter(SiteID=="Mosque")
#Create mosque only layer
grouped_building_plot <- grouped_building_plot %>%
addCircleMarkers(lng=mosque_only_dataframe$coords.x1,
lat=mosque_only_dataframe$coords.x2,
radius = 2,
stroke=FALSE,
fillOpacity= 1,
color = "#f19452",
label=~destruction_type,
group="Mosque") %>%
addLayersControl(overlayGroups = c("General Building", "Field", "Road", "School/University", "Industrial Facility", "Mosque"))
# Dataframe with only government buildings
government_only_dataframe <- coalition_strikes_2017 %>%
filter(SiteID=="Government Building")
#Create government building only layer
grouped_building_plot <- grouped_building_plot %>%
addCircleMarkers(lng=government_only_dataframe$coords.x1,
lat=government_only_dataframe$coords.x2,
radius = 2,
stroke=FALSE,
fillOpacity= 1,
color = "#e66a4d",
label=~destruction_type,
group="Government Building") %>%
addLayersControl(overlayGroups = c("General Building", "Field", "Road", "School/University", "Industrial Facility", "Mosque", "Government Building"))
# Dataframe with only hospitals
hospital_only_dataframe <- coalition_strikes_2017 %>%
filter(SiteID=="Hospital")
#Create hospital only layer
grouped_building_plot <- grouped_building_plot %>%
addCircleMarkers(lng=hospital_only_dataframe$coords.x1,
lat=hospital_only_dataframe$coords.x2,
radius = 2,
stroke=FALSE,
fillOpacity= 1,
color = "#d43d51",
label=~destruction_type,
group="Hospital") %>%
addLayersControl(overlayGroups = c("General Building", "Field", "Road", "School/University", "Industrial Facility", "Mosque", "Government Building","Hospital")) %>%
addLegend(pal = building_palette,
values = c("General Building", "Field", "Road", "School/University", "Industrial Facility", "Mosque", "Government Building","Hospital"),
title = "Type of building",
position = "bottomright") %>%
leaflet.extras::addResetMapButton()As the visualisation demonstrates the vast majority of buildings (10,732) are simply labelled “Building (General / Default)” and I wondered whether it was possible to find out further details using reverse geocoding. Using the GeoCage API, I looked up the latitude and longitude coordinates for each damaged building and extracted the formatted component of the response object which in some cases contained additional address details. After I exhausted my daily request limit (2500 queries) I noticed there were only 99 unique entries which didn’t contain the unhelpful “unnamed road” text so I didn’t add this new information to the dataframe. However it might to be interesting to explore whether other APIS could provide additional data. and if there are other fields of the GeoCage API which could be useful.
library(httr)
library(jsonlite)
library(dplyr)
### --- Load Raqqa Geographic data and select columns of interest
#Load UNOSAT data in dataframe form
raqqa_data <- readRDS("Datasets/Reshaped UNOSAT data.rds")
raqqa_data_summarised <- raqqa_data %>%
select(SiteID,Neigh,damage_year,destruction_type,coords.x1,coords.x2)
###--- Connect to GeoCage API perforl reverse geocoding on all destruction sites to get
# more information on the sites as just referred to as Building "(General / Default)" in the UNOSAT data
#Convert latitude and longitude columns to vectors
longitude <- (raqqa_data$coords.x1)
latitude <- (raqqa_data$coords.x2)
# Connect to OpenCage API by defining API key and endpoint
api_key <- readRDS("Secrets/api_key_open_cage.rds")
base_url <- "https://api.opencagedata.com/geocode/v1/json"
list_extra_info <- list() # for storing list of additional address information
for(x in 1:length(longitude)){
Sys.sleep=0.1 #to respect API limits
current_location <- paste(latitude[[x]],"+",longitude[[x]], sep="") # create query string from dataframe information
geocoding_response <- httr::GET(url = base_url,
query = list(key = api_key, q = current_location))
current_address <- httr::content(geocoding_response)
current_extra_info <- current_address$results[[1]]$formatted # additional information in formatted section of results
list_extra_info <-append(list_extra_info,current_extra_info)
}
#Convert list back to a dataframe
extra_data_dataframe <- data.frame(list_extra_info)
transpose <- data.frame(t(extra_data_dataframe))
rownames(transpose) <- seq(from=1, to=length(list_extra_info), by=1)
#Filter for useful additional information by dropping all rows that contain only "unnamed road"
filtered_address_info <- transpose %>%
filter(!str_detect(t.extra_data_dataframe., "unnamed road")) %>%
distinct()
#Export final list of additional information
saveRDS(filtered_address_info,"Visualisations/filtered address info.rds")library(rio)
library(data.table)
address_info <- rio::import("../Visualisations/filtered address info.rds")
address_info_table <- data.table(head(address_info, 5))
address_info_table## t.extra_data_dataframe.
## 1: Forosiyah Street, Raqqah Subdistrict, Raqqa, Syria
## 2: Mansour Street, Salhiyeh, Raqqa, Syria
## 3: فرن تشرين, 6, Baath, Raqqa, Syria
## 4: Makef Street, Batani, Raqqa, Syria
## 5: New Bridge Highway, Hettin, Raqqa, SyriaMedium Resolution Satellite Data (Sentinel-1 and Sentinel-2)
Overview
Launched by the European Space Agency in June 2015, primarily for agricultural monitoring and emergency management, Sentinel-2 collects land and ocean imagery to a maximum resolution of 10m every 5 days. The imagery consists of 13 spectral bands ranging from visible (red, green, blue) to near infrared (NIR) and short wave infrared (SWIR), it also includes specific frequencies for identifying vegetation (bands 5-8a) and clouds (Band 10: Cirrus).(“Custom Scripts Tutorial,” n.d.)
Sentinel -1 was launched by the European Commission’s Copernicus Project a year earlier in April 2014, and collects Synthetic Aperture Radar (SAR) data every 10 days. While Sentinel 2’s optical sensors allow it to capture imagery from reflected sunlight, Sentinel-1 produces its own electromagnetic beam (with 8 bands between 0.3GHz and 40GHz) and records its reflection. The longer wavelength radar frequencies not only pass straight through clouds so can be used in any weather, but penetrate soil and sand revealing hidden underground artefacts. Additionally, phase data can be extracted from the radar signal and used to calculate differences in surface topography with an accuracy of up to a few centimetres. Combining these height difference calculations with machine learning techniques, one 2022 studyidentified buildings damaged by earthquakes.
Real colour images can be constructed by superimposing red, green and blue bands but often it is more revealing to use false colour. The optimal frequency depends on the exact nature of the task. Infrared is generally preferable for forestry mapping as soil and vegetation are difficult to differentiate in true colour imagery (as their spectral reflectances are almost identical in the visible range), however there is a difference of almost 30% in the mid-infrared band. (“Custom Scripts Tutorial,” n.d.)
Frequencies can also be combined to accentuate features of interest, blending green (band 3), red (band 4) and infrared (band 8) creates a “false colour” image amplifying the infrared band. In 2020, Amnesty International used false colour imagery to identify vegetation burnt by gun propellant to pinpoint the firing position of Tigray strikes on Eritrea. Hundreds of standard combinations exist with the majority used to highlight vegetation, its water content and the chemical composition of the surrounding soil for use in crop monitoring.The images below show the bombing of Raqqa on October 10th, 2017 in different frequency bands. The real colour image clearly shows the plume of smoke over the city, but the SWIR band (based on bands 12, 8A and 4) shows individual fires. The Normalised Difference Moisture Index (NDMI) appears bright red over the smoke indicating the low level of humidity here, the bluer sections at the edges are fields so have a higher level of moisture.
Images from October 10 2017 in different frequency bands, showing real colour (left), SWIR (centre) and moisture index (right)
To choose suitable bands for identifying building destruction, I followed a simplified version of the methodology described by Pfeiffle (2022) and Putri et al. (2022). (Sandhini Putri, Widyatmanti, and Umarhadi 2022) [Conflict identification] Pfeiffle used machine learning methods to automatically identify building damage from airstrikes in Iraq and Northern Syria, choosing all Sentinel-2 images with a 10m resolution (red, green, blue and NIR). Assuming the damage profile of bombed buildings would be similar to those damaged in earthquakes, I also downloaded Sentinel-1 GRD imagery as Putri et al. concluded a fusion of Sentinel-1 and Sentinel-2 data yielded the highest accuracy in their random forest model identifying building damage from the 2018 Lombok earthquake. It would also be interesting to explore the possibility of using phase or polarisation parameters in the Sentinel-1 data to pinpoint destruction.
Data Download Options
Sentinel Satellite data is available free online through the following portals:
The EO browser is an interface allowing registered users to easily select and download satellite imagery of regions of interest from an interactive map. Using simple drop down menus users can also select the relevant time range, level of maximum cloud coverage and type of image. Results are presented as a simple image list which not only offers the option of fusing single images into standard combinations including real colour, false colour and Normalised Difference Vegetation Index (NDVI), but also of creating custom combinations. Images can be saved, shared, pinned, downloaded and automatically woven into time lapse sequences.
Alternatively Sentinel-1, 2 and 3 imagery can be downloaded via the Copernicus Hub API. The R package sen2r has been written to simplify the process of accessing Sentinel-2 data, allowing users to access imahergy not only via the libraries custom functions but also through a Shiny backed GUI interface. The RESTful Sentinel Hub Process API provides full access to data from all Sentinel satellites including polarisation and phase information. Users are offered the ability to both select products and perform bespoke image processing using the platform’s Javascript based Evalscript. (Example of commented evalscript from code)
Although it’s possible to easily visualise, download and customise satellite imagery using EO Browser and sen2r, this project downloaded all imagery directly from the Sentinel Hub API. To use satellite data in a potential machine learning project I wanted to create an efficient data pipeline to collect and process large numbers of satellite photos without scrolling through multiple images manually via the browser. The sen2R package does offer this functionality but only for Sentinel-2 data, also it requires sci hub log-in which is only activated one week after registration (and would have made meeting this deadline impossible :))
Accessing the Sentinel Hub API
Configuring the OAUTH 2.0 authorisation in R was a challenging task,
first it was necessary to register for a client ID and secret then
create an oauth app to generate a token (which is valid for 1 hour) The
Sentinal API requires a POST request specifying the endpoint and request
body. The{r, eval=FALSE} add_headers argument must also be
set, requiring: png images are returned
{r, eval=FALSE} Accept;the body request will be sent in
JSON {r, eval=FALSE} Content-Type and specifying the
necessary authorisation parameters.
library(httr)
library(brio)# for loading text file with no extension
library(rio)
library(png)
### --- Oauth Authorisation to access API
# IDs and secrets
client_id <- readRDS("Secrets/oauth_client_id.rds")
client_secret <- readRDS("Secrets/oauth_secret.rds")
# Create the app
app <- oauth_app(appname = "SatelliteImageryAnalysis",
key = client_id,
secret = client_secret)
# Create the oauth endpoint
endpoint <- oauth_endpoint(
request = NULL,
authorize = NULL,
access = "token",
base_url = "https://services.sentinel-hub.com/oauth")
# Cache the token to prevent API being called multiple times (token lasts 1 hour)
options(httr_oauth_cache=T) # prevents pop-up asking to cache token
token <- oauth2.0_token(endpoint = endpoint,
app = app,
use_basic_auth = T,
client_credentials = T)Using a combination of these parameters coupled with these examples it was possible to structure the following request in JSON format. The upper portion of the JSON specifies the date range, region of interest (bbox); acceptable level of cloud cover (maxCloudCoverage);type of satellite footage required (type), resolution (width, height) and coordinate system (crs) while the lower section consists of evalscript to process the imagery.
request={
"input": {
"bounds": {
"properties": {
"crs": "http://www.opengis.net/def/crs/OGC/1.3/CRS84"
},
"bbox": [
13.822174072265625,
45.85080395917834,
14.55963134765625,
46.29191774991382
]
},
"data": [
{
"type": "sentinel-2-l2a",
"dataFilter": {
"timeRange": {
"from": "2018-10-01T00:00:00Z",
"to": "2018-12-31T00:00:00Z"
}
}
}
]
},
"output": {
"width": 512,
"height": 512
}
}Ideally the POST request should be sent as a multipart form, which would allow the evalscript to be written in ordinary JSON. HoweverI was unable to do this correctly, so as a workaround escaped the entire JSON Evalscript section (using this tool) and sent it as part of the same request. The image information was then extracted from the response object and the file saved as png.
### --- Build request
request <- read_file("R Scripts/POST body request") # loads request from JSON in text file
sentinel_endpoint <- "https://services.sentinel-hub.com/api/v1/process"
# Specify which formats are accepted in request(JSON) and response(png)
response <- httr::POST(sentinel_endpoint,
body=request,
add_headers(Accept = "image/png",
`Content-Type`="application/json",
Authorization = paste("Bearer", token$credentials[[1]], sep = " ")))Function to automatically download Sentinel-2 Data
To simplify the download process I converted this code to a function to download the Sentinel-2 data automatically. The function handles the authentication process and accepts the following arguments:
download_satellite_data(start_date, period_length_days, number_periods, min_lat,max_lat,min_long,max_long,satellite_type,cloud_cover,pixels_dimensions_height)
# Minimum and maximum latitudes and longitudes specify region of interest
# period_length, start_date and number of periods can be combined to specify data collection period
# For example: download_satellite_data("2017-05-01", 7, 32)
# Would return the best cloud free image taken every week from May 2017- September 2017 (32 in total)
# pixels_dimensions_height specifies the horizontal dimension of the output image in pixelsTo make the file structure as clear as possible the JSON body of the request was stored in the separate text file (“POST body request.txt”) and imported to the R file containing the function (“Download Sentinel Satellite Data.R”. Str_replace() was used to replace the existing parameters with the user defined function arguments in the JSON code. str().
request_with_start <- str_replace(request, "InsertStartDateHere", start_date_as_string)
request_with_end <- str_replace(request_with_start, "InsertEndDateHere", end_date_as_string)
### Create satellite type parameter of request body----------
request_with_satellite <- str_replace(request_with_end, "InsertSatelliteTypeHere", satellite_type)
### Create cloud cover parameter of request body----------
request_with_cloudcover <- str_replace(request_with_satellite, "InsertCloudCoverHere", as.character(cloud_cover))It is easiest to define the area of interest using the maximum and minimum longitudes and latitudes, however using the bbox argument in the JSON selected a much larger area than expected, so this function uses the type Polygon parameter instead. This argument takes a series of geo-coordinates and joins them up so the first and last coordinate must be the same to close the shape.
It was also necessary to specify the dimensions of the returned images as selecting width and height at random led to distorted, skewed pictures. Using the code below the size of region of interest in degrees was calculated then converted to kilometres using standard formulae. By calculating the ratio of the sides of the rectangle, the optimum height width ratio was calculated and combined with the user’s choice of pixels_dimensions_height to define the optimal image dimensions.
### Calculate optimum image size----
# API demands user specifies the dimensions of the output image
# To ensure it is not stretched, dimensions have been calculated based on user specifying desired number of pixels in horizontal direction
# Calculate longitude and latitude differences from user input
latitude_diff <- max_lat-min_lat
longitude_diff <- max_long-min_long
#Convert difference in latitude from degrees to km using standard conversion formula
lat_metres <- abs(111.32*latitude_diff)
long_metres <- abs(longitude_diff*40075*cos(max_lat)/360)
# Calculate ratio of height:length of the selected satellite image
normalised_lat_km <- lat_metres/long_metres
normalised_lat_km
# Calculate width of output image based on user defined pixel height
pixels_dimensions_width= round(normalised_lat_km*pixels_dimensions_height)Finally the date the image was taken was added to the image, this was simple as the date (current_date) was already saved in memory. The function was then used to download 32 real colour satellite images from May-November 2017 with a maximum cloud cover of 10%.
### Annotate png with Date---
#Extract date parameters to display
year <- lubridate::year(current_date_start)
month <- lubridate::month(current_date_start, label=TRUE)
day <- lubridate::day(current_date_start)
complete_date <- paste(day, month, year, sep=" ")
processed_png <- image_read(data_filename)
annotated_png <- image_annotate(processed_png, complete_date,location = "+10+450", size = 30, color = "white")Data Visualisation: Coalition Airstrikes Against Raqqa: May-October 2017
Using this tutorial an animated GIF was created from the downloaded satellite imagery showing a stop animation of the 4-month bombardment of Raqqa.The animation clearly shows the destruction of the city starting in the suburbs in the North West of the city and gradually moving towards the centre. On October 9th a huge plume of black smoke is clearly visible.
